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Abstract 

We present a numerical method for scalar conservation laws in one 
space dimension. The solution is approximated by local similarity so- 
lutions. While many commonly used approaches are based on shocks, 
the presented method uses rarefaction and compression waves. The so- 
lution is represented by particles that carry function values and move 
according to the method of characteristics. Between two neighboring par- 
ticles, an interpolation is defined by an analytical similarity solution of 
the conservation law. An interaction of particles represents a collision 
of characteristics. The resulting shock is resolved by merging particles 
so that the total area under the function is conserved. The method is 
variation diminishing, nevertheless, it has no numerical dissipation away 
from shocks. Although shocks are not explicitly tracked, they can be 
located accurately. We present numerical examples, and outline specific 
applications and extensions of the approach. 

1 Introduction 

Hyperbolic conservation laws are important models for the evolution of con- 
tinuum quantities, and their efficient numerical approximation still involves 
challenges. A fundamental property is the occurrence of shocks. These ap- 
pear when parts of the solution that move at different characteristic velocities 
collide, forming a traveling discontinuity. A sharp and accurate resolution of 
shocks, without creating spurious oscillations, is one crucial challenge in nu- 
merical methods. Traditional finite difference [T7] or finite volume [5] methods 
operate on a fixed grid. Straightforward linear numerical schemes are either 
of low accuracy, or create overshoots near shocks [5] . These overshoots can be 
avoided by considering more complex nonlinear schemes, such as finite volume 
methods with limiters [24], and ENO [9] or WENO [20] schemes. These difficul- 
ties show up in simple problems, such as scalar, one dimensional conservation 
laws, even for linear advection. One interpretation of this challenge is that fixed 



grid methods attempt to represent advection by a change in function values. 
However, since advection moves the solution along the a;— axis, it is natural to 
move the position of the computational nodes instead. Two popular approaches 
that employ this either track shocks explicitly, or trace characteristics. One of 
the first shock-tracking methods is Godunov's scheme [5]. An initially piecewise 
constant function is evolved analytically by solving the local Riemann problems 
exactly. This works well until two shocks (or rarefactions) interact. Godunov's 
scheme avoids the resolution of this interaction by remeshing: the solution is 
averaged over cells, and a new piecewise constant solution is reconstructed. In 
this sense, Godunov's scheme is actually a fixed-grid finite volume method, and 
thus incurs the problems described above. The idea of tracing characteristics 
has been formulated by Courant, Issacson and Rees in the CIR method [3]. The 
solution of a conservation law is given by a simple ordinary differential equation 
along a characteristic curve. The CIR method updates function values at grid 
points by following characteristic curves. Function values between grid points 
are defined by an interpolation, which is a remeshing step. Therefore, the CIR 
method is a fixed grid approach as well. In the simplest case of using linear 
interpolation between grid points, a simple non-conservative upwind scheme is 
obtained. 

Most popular numerical schemes for conservation laws are fixed grid ap- 
proaches. Advantages include simple data structures and an easy generalization 
to higher space dimensions using dimensional splitting. While these methods 
can be derived from similarity solutions (finite volumes) or the method of char- 
acteristics (CIR), they require a global remeshing step, which introduces numer- 
ical errors, numerical dissipation and dispersion, and a decay in entropy even 
away from shocks. More recent approaches attempt to avoid global remeshing, 
or at least weaken its negative impact. An example of such an approach is 
front tracking [13] . Shocks are tracked explicitly, but unlike with Godunov-type 
schemes, no cells exist and no averaging is performed. Instead, the interaction 
of the local similarity solutions is resolved where necessary. The difficulty of re- 
solving interactions of rarefactions with shocks is avoided by approximating the 
flux function by a piecewise linear function. Thus, only discontinuities have to 
be tracked, and rarefactions are approximated by step functions. Analogously, 
the method of characteristics can be used wherever characteristic curves do not 
intersect. Lagrangian particle methods are based on this idea. Particles follow 
the characteristics, and only when they come too close is their interaction re- 
solved locally. Classical approaches, such as smoothed particle hydrodynamics 
[2"2"] . consider a pressure that prevents particle collisions. While such approaches 
work well in smooth regions, a stable treatment of shocks is challenging. In ad- 
dition, hybrid methods exist, such as the self-adjusting grid method |l(Jj . which 
operates on a fixed grid, and, in addition, tracks up to one shock particle in 
each cell. 

The method presented in the current paper can be interpreted as a combi- 
nation of some of the aforementioned ideas. Characteristic curves are traced 
using Lagrangian particles. In addition, between neighboring particles, a sim- 
ilarity solution is evolved exactly. While front-tracking methods trace shocks, 
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the current method uses rarefaction and compression waves to represent the 
solution. Due to this analogy, the approach can be called rarefaction tracking. 
The local similarity solutions give rise to a natural definition of area between 
two particles. This, in turn, admits a merging strategy of colliding particles 
that conserves area exactly. Thus, shocks are represented by compression waves 
that are removed upon breaking. Due to the correct evolution of area, correct 
shock speeds are obtained. The method was originally motivated in the context 
of Lagrangian particle methods [5J. Well-posedness, accuracy and convergence 
have been shown and investigated in [7J. This paper focuses on the connection 
to local similarity solutions. In addition, applications for which the approach 
has advantages over fixed-grid approaches are outlined. The method, as it is 
presented here, is suited for scalar conservation laws, since for this kind of equa- 
tions similarity solutions are straightforward to construct. Also, only one space 
dimension is currently considered, since in this case a shock can only happen 
if particles collide. However, extensions to balance laws with source terms are 
presented, and other possible generalizations are outlined. 

In Sect. [2J we outline the general characteristic particle approach. In Sect. [31 
similarity solutions between particles are constructed. These give rise to an ex- 
act expression for area, which is the basis for particle management, as explained 
in Sect. [U Properties of the approach as presented in Sect. and fundamen- 
tal extensions (for simple sources, and flux functions with one inflection point) 
are outlined in Sect. [BJ In Sect. [71 results of computations with the presented 
method are given. We consider an "academic" flux function, and two phys- 
ical models. Possible applications of the method and possibilities for further 
extensions are presented in Sect. [U 



2 The Particle Method 

We consider a one dimensional scalar conservation law 

u t + f{u) x = 0, u(x,0)=u o {x) (1) 

with /' continuous. In addition, until we explicitly relax the condition, we as- 
sume that on the entire interval of considered function values [min^. u(x), max x u(x)], 
/ is cither convex or concave. The characteristic equations [I] 

U = f'(u) 

\u = y ' 

yield the solution (while it is smooth) forward in time: at each point (xq, uq(xo)) 
a characteristic curve x(t) = xq + f'(uo(xo))t starts, carrying the function 
value u(x(t),t) = u (x ). The method of characteristics gives rise to a particle 
method: 

1. Sample the initial function uo(x) by a finite number of points (xi, u,). The 
aspect of sampling the initial conditions is addressed in Sect. 13.21 
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2. The solution at time t > is given by the points (xi + f'(ui)t, m). At these 
points, the solution is represented exactly, without any approximation 
error. 

This basic approach incurs two fundamental limitations, which need to be reme- 
died before it can become a useful method: 

• The solution is only known at the particle positions. If neighboring par- 
ticles depart, large regions of unspecified function value may arise. The 
presented approach remedies this problem by defining an interpolation be- 
tween particles. This allows an easy insertion of new particles into gaps 
wherever desired. 

• A particle may overtake its neighbor, corresponding to the intersection of 
characteristic curves. The correct weak solution of the conservation law JT]) 
develops a shock, i.e. a moving discontinuity, into which the characteristic 
curves keep running [4|. In the presented approach, particles are merged 
upon collision. An interpolation between particles allows the merging of 
particles in such a way that area is conserved. 

The insertion and merging of particles is called particle management. It is 
described in detail in Sect. 01 Particle management is the key ingredient that 
enables us to use the method of characteristics beyond the first occurrence of 
shocks. 

Given the solution at time t, represented by particles (Xi(t), Ui(t)), the time 
of the next particle collision is t + At s , where 

At- - Xi+1 ~ Xi 

l ~ f'(u i+1 )-f'( Ui ) [6) 

is the time until the collision between particles i and i + 1, and 

At s = min{{Ati\Ati > 0} Uoo} (4) 

is the time until the next collision in the future. We can step forward by up 
to At s without risking a collision. After At s , at least one particle will share 
its position with another. To proceed further, we need to perform particle 
management. We insert new particles into gaps between particles of distance 
larger than a maximum distance d max - Then, each pair of particles sharing their 
position is merged. Optionally, one can perform a merging step for all particles 
closer than a minimum distance d m ; n . Choices for the maximum and minimum 
distance are described in Sect. |U 



3 Interpolation by Similarity Waves 

Consider two neighboring particles (xi(t),ui), (x2(t),U2), with X\ < Xi- If 
/'(tti) > fiu-i), the two particles meet in the future, at time t + At\, where 
At\ > is given by ([3]). Let us assume that between the two particles, the 
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solution is continuous at all times before t + At\. This implies that the meeting 
of the particles is the first occurrence of a shock, and at this time the solu- 
tion is a discontinuity at position x^ = x% + f'(ui)Ati. Any point (x s h,u), 
u G [1*1,112] on the shock can be traced backwards along its characteristic to 
(xgh — f'(u)Ati,u). As shown in Fig.[Tl this defines the interpolation uniquely 
as 

x(u) =xx + f'[ui)Atx - f'(u)At 1 =xt+ {} U \ — Krr—s {x 2 - an) . (5) 

f'{U2) - /'(Ml) 

In the case of departing particles, i.e. f'{u\) < f'(u 2 ), the above construction 
is exactly reversed in space and time. We assume that the solution between 
the particles came from a discontinuity at some time in the past. Consequently, 
formula §5§ also describes this case. 

Theorem 1. The interpolation ([5]) between two particles is a solution of the 
conservation law ([T]) . Consequently, for multiple particles, the function defined 
piecewise by the interpolation ^ is a classical (i.e. continuous) solution to the 
conservation law (TTJ), until particles collide. 

Proof. Using that Xi(t) = f'(ui) for i = 1, 2 (from {2j) one obtains 

dx . f(u)-f( Ul ) . 

di^ t) = Xl+ f'(u 2 )-f'(u^ X2 - Xl) 

= fM + / 11) - (//(M2) - fM) = f{u) • 

Thus every point on the interpolation u(x, t) satisfies the characteristic equation 
((2]). At the particles, the solution possesses kinks. These can be interpreted as 
shocks of height zero, and thus satisfy the Rankine-Hugoniot conditions [4]. 
Hence, the full piecewise wave solution is a weak solution of |TJ . □ 

The interpolation ([5]) defines either a rarefaction wave (if f'(ui) < /'(^a)), 
or a compression wave (if f'(ui) > /'(M2)) between two particles. In the case 
Mi =1*2, a constant interpolant is defined. Hence, at all times the solution is 
approximated by similarity waves. 

Remark 1. That © defines the interpolant as a function x(u) is, in fact, an 
advantage, since at a discontinuity x% = x%, characteristics for all intermediate 
values u are defined. Thus, rarefaction fans arise naturally from discontinuities 
if f'{u\) < f'iu-i). If / has no inflection points between u\ and 1x2, the inverse 
function u(x) exists. For plotting purposes we plot x(u) instead of inverting the 
function. 

3.1 Evolution of Area 

The interpolation defines an area between particles. As shown in [7], the inte- 
gration of ([5|) yields 

x 2 (t) 

u(x,t)dx = (x 2 (t) - x 1 {t))a{u 1 ,u 2 ) , (6) 

xi(t) 
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where a(u\, u 2 ) is the nonlinear average function 



[f(u)u f{u)]Z Cf"(u)udu 
[J W\ U1 J Ul f"(u)du 

The integral form shows that a is indeed an average of u, weighted by /". In 
[7 it is shown that 0(1*1,1*2) is symmetric, bounded by its arguments, strictly 
increasing in both arguments, and continuous at u\ = u 2 . 

The area under the interpolant can also be derived from the conservation 
law by observing that the change of area between two characteristic particles 
(xi(t),ui) and (x 2 (t),u 2 ) is given by 

d f X2< -^ 

— u(x,t)dx=(f'(u 2 )u 2 -f(u 2 ))-(f'(u 1 )u 1 -f(u 1 )) 
dtJ Xl{t ) (8) 

= F{u 2 )-F{u 1 ) = [F\l\ , 

where F = f (u)u — / is the Legendre transform of the flux function /. Since 
u is conserved along the characteristics, F is conserved as well. Thus, the area 
between particles changes linearly in time, as does their distance 

j t {x 2 {t) - Xl {t)) = x 2 {t) - xt(t) = f\u 2 ) f'( Ul ) = [f'{u)} u u \ . (9) 

As before, these results require that the solution remains continuous until the 
particles collide. At the time of collision to, both the distance and the area 
vanish. Thus ([Sj) and ([9|) integrate to 

x 2 (t) 

u(x,t)dx = (t-t )-[F(u)]2 , and 

x x (t) 

X2 (t)-x 1 (t) = (t-t ).[f(u)C i , 

which lead to ©. 

Remark 2. The definition of area between particles gives rise to an exactly 
conservative resolution of particle interaction. This poses an interesting analogy 
to finite volume methods. While these are derived by the flux / through fixed 
boundaries 

d r u (x,t)dx =-mz , 



dt 



the presented method is based on the Lagrangian flux F through moving bound- 
aries, as given by (8J. 

3.2 Sampling of the Initial Conditions 

Given particle positions and function values, the interpolation ((5|) defines a func- 
tion everywhere. If the initial conditions uq(x) can be represented exactly by 
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Figure 1: Definition of the interpola- 
tion. After time Ati, particles 1&2 
form a shock. From this x{u) is found. 



Figure 2: Merging particles 2&3. The 
values (x23,U23) are found so that the 
area of the two "triangles" is equal. 



finitely many particles and the interpolation ([5]) in between, then the charac- 
teristic evolution yields an exact solution of the conservation law. In general, 
the initial conditions must be approximated by the piecewise interpolation ([5]). 
This approximation step is similar to finite element approaches. The selection 
of initial particle positions Xi corresponds to the construction of a mesh. The 
function values Ui can be chosen to minimize the approximation error between 
the exact initial function and the numerical approximation. However, there are 
two aspects of the initial sampling that are very specific to conservation laws 
and the numerical approach: 

• The initial function u (x) may involve discontinuities. We assume only 
finitely many jumps, and that they are known. Each jump can be repre- 
sented exactly by two particles at the same position (but with different 
function values). Having this option is one definite advantage of particle 
methods over Eulerian approaches. 

• Area is of great importance for conservation laws. The interpolation yields 
an exact knowledge of area in the numerical approximation, and the ap- 
proximation step of the initial function can incorporate this aspect, and 
chose an initial sampling that respects the area of uq(x). 

A simple approach is to represent jumps exactly, and then sample positions Xi 
equidistantly between jumps while choosing function values Uj = Uo(xi). Ob- 
viously, better sampling strategies are possible, but this simple approach often 
yields satisfying results. In [7], we prove that with this simple sampling strategy, 
the piecewise interpolation (|5|) is a second order accurate approximation to any 
C 2 function with \ f" (uq(x))\ > C > 0, i.e. the initial function does not cross 
nor touch any inflection point of /. For non-convex flux functions, as outlined 
in Sect. 16.11 the second order accuracy is, in general, lost when uq crosses an 
inflection point of /. 
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4 Particle Management 



The time evolution of the conservation law (fTJ) is described by the characteristic 
movement of the particles ([2]) . Particle management is an "instantaneous" op- 
eration (i.e. happening at constant time) that allows the method to continue an 
evolution forward in time. It is designed to conserve area: The function value 
of an inserted or merged particle is chosen so that the total area is unchanged 
by the operation. 

Consider four neighboring particles located at x\ < xi < X3 < xJ3 with as- 
sociated function values u±, u 2 , 113, 114, and u 2 =£113. There are two possible 
cases that require particle management between Xi and X3: 

• Insertion: If two particles deviate from each other (f'(u 2 ) < /'(U3)), 
their distance may become unacceptably large: X3 — x-i > d max . In this 
case, we insert a new particle (223,^23) with X2 < £23 < X3, such that the 
area is preserved: 

(x 2 3 - x 2 ) a(u 2 , U23) + (x 3 - x 23 ) a("23, ""3) = (x 3 - x 2 ) a(u 2 , u 3 ) . (10) 

One can, for example, set £23 = X2 + Xa and find U23 by (fTU)) . or set u 2 3 = 
and find X23 by ([10]). 

• Merging: If two particles approach each other (f'(u 2 ) > /'(U3)), they 
will eventually collide: X3 — x 2 < d m i n . In this case, we replace both with 
a new particle (£23,1*23) a ^ x 23 ~ X2 + X3 , with M23 chosen, such that the 
area is preserved: 

(X23 - Xi) a(ui, U23) + (X4 - X23) a(u23,Ui) 

= (x 2 -Xi)a(ui,u 2 ) + {X3-X2) a(u 2 ,u 3 ) + (x 4 -X3) a(u 3 , u 4 ) . 

Figure [2] illustrates this merging step. The merging step is extended by 
an entropy fix, as described in Sect. 14. H 

The two distance parameters have very natural meanings for the method. The 
maximum distance d max is the resolution of the method. Due to insertion, the 
distance between particles on rarefactions is always between d max /2 and c? max . 
The minimum distance <i m i n is ideally zero, therefore it is not a real parameter 
of the method. For stability reasons, one can assign a small positive value. 
This is of interest in extensions to balance laws, in which case the characteristic 
equations may have to be approximated numerically. 

We observe that inserting and merging are similar in nature. Conditions (1101) 
and (jlip for M23 are nonlinear (unless / is quadratic, see Rem. |3]). For most 
cases, U23 = " 2 ^" 3 is a good initial guess, and the correct value can be obtained 
(up to the desired precision) by a few Newton iteration steps (or bisection, if 
the Newton iteration fails to converge). It is shown in [7], that the new particle 

1 If more than two particles are at one position (x), all but the one with the smallest value 
(it) and the one with the largest value (it) are removed immediately. 
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value U23 for intersection and merging always exists, is unique, and bounded by 
the surrounding particle values 1*23 G [^2,^3]- For merging, the latter result is 
easy to show when X2 = X3, but it also holds when the slope \us — 1*2 1/ 1^3 — £2 
is larger than a given threshold (which depends on the surrounding particle 
geometry x\, . . . , X4, u\, . . . , U4). Thus, the merging step is robust with respect 
to small deviations in the distance of the merged particles. This also holds for 
the case xi > X3, i.e. merging can in principle be performed after two particles 
have overtaken. 

4.1 Entropy Fix for Merging 

The merging of particles resolves shocks, which are moving discontinuities. A 
weak formulation of the conservation law |T]) admits these solutions. For some 
initial conditions, the weak solution concept allows for infinitely many solutions, 
and an additional condition has to be imposed to single out a unique solution. 
This can be done by defining a convex entropy function r\, and an entropy flux 
q, and requiring the entropy condition 

7](u) t + q(u) x < (12) 

to be satisfied. This condition is equivalent to obtaining shocks only when 
characteristics run into them [2J. It is shown in [TSJ Chap. 2.1] that if (|12p 
is satisfied for the Kruzkov entropy pair r)k(u) = \u — fc|, qk(u) = sign(w — 
k)(f(u) — f(k)), for all k, then it is satisfied for any convex entropy function. 

Relation (fT2l) implies that the total entropy / r](u(x, t)) dx of a solution 
does not increase in time. Particle merging changes the solution locally. It is 
designed to conserved area exactly. Similarly, we require that the merging step 
not increase the total entropy. In [7], it is shown that this property is satisfied 
if shocks are reasonably well resolved. 

Lemma 1. Let x\ < xi — X3 < x^ be the locations of four particles, with 
particles 2 and 3 to be merged, thus ui > U3 ( if f" > 0), respectively ui < 1*3 
(if f" < 0). If the value U23 resulting from the merge satisfies u± > 1/23 > Ui 
(if f" > 0/, respectively U\ < U23 < (if f" < 0), then the Kruzkov entropy 
J \u — k\ dx does not increase due to the merge. 

This entropy condition is ensured by an entropy fix: A merge is rejected 
a posteriori if the condition is not satisfied. In this case, points are inserted 
half-way between the shock and the neighboring particles, and the merge is re- 
attempted. The condition is satisfied if new particles are inserted sufficiently 
close to the shock, hence this approach is guaranteed to succeed. 

4.2 Shock Postprocessing 

The particle method does not explicitly track shocks. However, shocks can be lo- 
cated. Whenever particles are merged, the new particle can be marked as a shock 
particle. Thus, any shock stretches over three particles [x\, ui), (x2, U2), (^3, U3), 
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with the shock particle in the middle. Before plotting or interpreting the solu- 
tion, a postprocessing step can be performed: The shock particle is replaced by 
a discontinuity, represented by two particles (2I2, Ui), {%2, U3), with their posi- 
tion X2 chosen such that area is conserved exactly. This step is harmless, since 
an immediate particle merge would recover the original configuration. The nu- 
merical results in Sect. [73] indicate that, with this postprocessing, the resulting 
solution is second order accurate in L . In particular, the shock is located with 
second order accuracy. 

5 Properties of the Method 

The presented method consists of two ingredients: time evolution by character- 
istic particle movement, and particle management at stationary time. Due to 
the concept of similarity solutions between particles, an exact solution of the 
conservation law (fTJ) is evolved during particle movement (see Thm. [T|). 

Theorem 2 (TVD and entropy diminishing). Solutions obtained by the particle 
method are total variation diminishing and satisfy the entropy condition (|12[) . 
Hence, no spurious oscillations are created, and correct entropy solutions are 
approximated. 

Proof. Due to Thm. [T] the characteristic particle movement yields an exact clas- 
sical solution. Thus, both the total variation and the entropy of the solution are 
constant. Particle insertion simply refines the interpolation, without changing 
the solution. Particle merging is shown to yield a new value bounded by the 
surrounding particle values (Sect. 2]), thus it is TVD. The entropy fix fSect. ETTj) 
ensures the entropy condition is satisfied. □ 

Remark 3 (Quadratic flux function). The method is particularly efficient for 
quadratic flux functions. In this case the interpolation ([5]) between two points 
is a straight line, and the average ((Jj is the arithmetic mean a(m, U2) — Ul + 112 . 
Thus, the function values for particle management can be computed explicitly. 

Remark 4 (Computational Cost). An interesting aspect arises in the computa- 
tional cost of the method, when counting evaluations of the flux function / and 
its derivatives /', /". Particle movement does not involve any evaluations since 
the characteristic velocity of each particle f'(ui) does not change. Consider a 
solution on t € [0, 1] with a bounded number of shocks, to be approximated by 
0(n) particles. Every particle merge and insertion requires O(l) evaluations. 
After each time increment ((4|), 0(1) management steps are required. The total 
number of time increment steps is 0(n). Thus, the total cost is 0(n) evalua- 
tions, as opposed to 0(n 2 ) evaluations for Godunov-type methods. Note that 
this aspect is only apparent if evaluations of /', /" are expensive, since the total 
number of operations is still 0{n 2 ). 
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6 Fundamental Extensions and Generalizations 



The basic method, presented in the previous sections, is formulated for scalar 
one-dimensional hyperbolic conservation laws with a flux function that is convex 
or concave. In this section, we present how some of these restrictions can be 
relaxed. In Sect. 16.11 we describe the treatment of flux functions with one 
inflection point. Section l(T2l addresses the incorporation of simple source terms. 

6.1 Non-Convex Flux Functions 

The key idea in generalizing the method for flux functions that have inflection 
points is to carry inflection particles, such that between neighboring particles, 
the flux function is always either convex or concave. At inflection particles, the 
interpolant ([5]) has an infinite slope (since /" = at this point), which poses 
some difficulty for a high order approximation of initial functions that cross an 
inflection point [TJ. The characteristic particle movement, and the expressions 
for area and average ([TJ apply as before. However, the merging of particles when 
an inflection particle is involved requires a special treatment. The standard 
approach, as presented in Sect. 21 replaces two particles by one particle of a 
different function value. If an inflection particle is involved in a collision, the 
merging must be done in a different way so that an inflection particle remains. 

The key idea is to move the inflection particle's position, while preserving its 
function value. Depending on the configuration, one of three different merging 
strategies has to be performed. Figure [3] shows these merging strategies for 
the case of a single inflection point (we do not consider here the interaction of 
multiple inflection points), with /"' > (the inflection particle is the slowest), 
and particles 2 and 3 to be merged. The other cases are simple symmetries of 
this situation. Each of the following strategies is attempted if the previous one 
failed: 

1. Remove particle 2 and increase X3 such that area is preserved. Accept if 
X3 need not be increased beyond X4. 

2. Remove particle 2, set X3 = X4 and increase both such that area is pre- 
served. Accept if new value for X3 and X4 is smaller than x$ . 

3. Remove particle 4, set X3 = x$ and lower ui such that area is preserved. 

It is shown in [TJ that one of these three options must be accepted. Note that 
the final configuration may involve a new discontinuity {X3 — 24 or X3 = X5). 
However, this is not a shock, but rather a rarefaction, and the particles will move 
away from each other. Consequently, these particles should not be merged. The 
presented strategy guarantees that in each merging step one particle is removed. 

6.2 Simple Source Terms 

Source terms extend the conservation law Q by a nonzero right hand side 

u t + f(u) x = Su , (13) 
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Figure 3: Particle management around an inflection particle (/"(M3) = 0) 




where S is an operator on u. Here we consider the case where the source term is 
a simple function Su = g(x,u), as opposed to a differential or integral operator. 
An example of S being an integral operator is outlined in Sect. 18.31 

If the source is a simple function Su = g(x, u), equation p"3|) becomes a bal- 
ance law. In this case, the method of characteristics [4] applies. The presented 
particle method can be generalized as follows. During the particle movement, 
the characteristic equations 

fx = f'(u) , x 

• 7\ (14) 

are solved numerically, for instance by an explicit Runge-Kutta scheme. Merging 
takes place when two particles are too close. Particle management is performed 
according to (|10[) and (jlip . as derived for the source- free case. This is an 
approximation, since the interpolation (J5|) is not a solution of the balance law 
(|13|) . If the advection dominates over the source, we expect the approximation to 
be very accurate. The numerical results in Sect. [731 indicate that this approach 
can yield good results even when the source and the advection term are of 
comparable magnitude. The solution is found to have second order convergence, 
and accuracy that compares favorably to CLAWPACK [T] . 
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Figure 5: L 1 Error of the particle method for f(u) = ju 4 , before and after 
a shock. In the left plot the error is only due to the initial sampling. In the 
right plot the advantage of "postprocessing" is made apparent, as it converts 
the first-order accuracy of the "pure" method, into second-order. 



7 Numerical Examples 

The particle method is numerically investigated in various examples. In all 
cases, the reference solution is obtained by a high resolution CLAWPACK [l\ 
computation. A more rigorous exposition of the order of convergence and com- 
parisons with other methods are presented in [7]. 

In Sect. 17. li the evolution and the formation of shocks for smooth initial 
data under a convex flux function are presented. The reduction of order due 
to a shock (if no postprocessing is done) is shown, as is the effectiveness of 
the postprocessing in recovering the second-order accuracy. In Sect. 17. 2\ the 
Buckley-Leverett equation is considered, as an example of a non-convex flux 
function. In Sect. 17.31 Burgers' equation with a source is simulated. The source 
code and all presented examples can be found on the particleclaw web page 



7.1 The Formation of Shocks 

Figure S] shows the smooth initial function uq(x) = cxp(— a; 2 ) cos(7rx), and its 
time evolution under the flux function f(u) = \u A . Note the curved shape 
of the interpolation |JSJ|. Initially, we sample points on the function uq. At 
time t = 0.25, the solution is still smooth, thus the particles lie exactly on the 
solution. At time t = 8, shocks and rarefactions have occurred and interacted. 
Although the numerical solution uses only a few points, it represents the true 
solution well. 

The error for different numbers of particles is shown in Fig. [5] Before shocks 
have emerged, the solution is smooth and the error is only due to the initial 
sampling, and therefore second-order. After a shock has emerged, the situation 
is not so simple: the exact solution has a jump discontinuity, and we measure 
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Figure 6: The particle method solution for the Buckley-Leverett equation at 
three different times. The exact solution is given in gray, and the dots show the 
numerical approximation. 



the error according to the L 1 norm of the difference between the interpolation 
and the "exact" solution. Therefore, since we normally have just one particle 
near the shock, we will make an error proportional to the height x width of the 
shock. Since the width of the shock is generally proportional to the resolution 
of the method, we get first order accuracy. However, since we have the correct 
area, the postprocessing step, explained in Sect. 14.21 recovers the second order 
accuracy of this solution, as is evident from the figure. 

7.2 A Non-Convex Flux Function 

As an example of a non-convex flux function, we consider the Buckley-Leverett 
equation, which is a simple model for two-phase fluid flow in a porous medium 
(see LeVeque [T51 Chapter 16]). The flux function is 

f(u)=u 2 /(u 2 + ±(l-u) 2 ). (15) 

We consider piecewise constant initial data with a large downward jump that 
crosses the inflection point u* — 0.387, and a small upward jump: 

!1 x < 1 
1< x < 1.3 (16) 
0.3 1.3 < x. 

The large jump develops a shock at the bottom and a rarefaction at the top, 
while the small jump develops a pure rarefaction. Around t — 0.2, the two sim- 
ilarity solutions interact, thereby lowering the separation point between shock 
and rarefaction. The solution is shown in Fig. [6] The analysis presented in [7] 
shows that the convergence obtained here is less than second order, due to the 
error on the intervals that abut the inflection particle. 
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Figure 7: Numerical results for Burgers' equation with a source. We can see that 
as time progresses, a bump is formed, a shock/rarefaction structure emerges, 
the wave passes over the "bump" , and the two shocks merge into one. 

7.3 A Simple Source Term 

We consider Burgers' equation with a source 



It is a simple model for shallow water flow over a bottom profile b(x). As in [16] . 
we consider the domain x S [0, 10]. We include the source into the method of 
characteristics, as explained in Sect. 16.21 Figure [7] shows the numerical solution 
for the initial condition 



The particle method (dots) approximates the exact solution (grey line). A 
particular aspect in favor of the characteristic approach is the precise (up to the 
precision of the Runge-Kutta scheme) recovery of the function values after the 
obstacle. 

8 Further Applications and Extensions 

The presented numerical method employs the method of characteristics and re- 
solves the interaction of characteristic particles. This implies important advan- 
tages: the scheme is exact away from shocks, the provided interpolation yields a 
certain level of "subgrid" resolution, shocks are located very accurately, and the 
method is TVD. A crucial drawback is that the basic method is limited to scalar 
problems in one space dimension. However, this restriction is not as limiting as 
it may seem at first glance. Many complex problems that arise in applications 
consist of, or can be reduced to, scalar one-dimensional sub-problems. As we 
motivate in the following, the presented particle method is a good candidate for 
the numerical approximation of such sub-problems. 
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8.1 Flow on Networks 



In many models of flows on networks, the evolution on each edge is governed by 
a scalar one-dimensional conservation law, with appropriate coupling conditions 
on the network nodes. An important example is the simulation of vehicular traf- 
fic on road networks. Optimization of traffic flow is a topic of current research 
[12] , and frequently a key challenge in such simulations is simply the size of the 
network. For large networks, only a few number of unknowns can be attributed 
to each individual road segment. Yet, exact conservation properties and accu- 
rate shock location are required in order to track traffic jams. In addition, traffic 
densities must never become negative, hence TVD methods are desirable. The 
presented particle method satisfies all those requirements: the method is TVD, 
the solution is represented accurately with just a few particles, and shocks are 
located accurately. 

A simple traffic model that is commonly considered for network flows is 
the Lighthill-Whitham model [15]. On each road segment, the evolution of the 
vehicle density p{x,t) is given by the scalar conservation law 

Pt + Mp)L - . (18) 

Here v(p) is the speed at which vehicles at density p move. It is a decaying func- 



tion. A popular choice is v(p) 



1 - -e- 

Pma: 



where v max is the speed limit 



assumed on an empty road, and p max denotes the maximal density of vehicles 
on the road. This choice of v(p) yields a quadratic flux function, for which the 
presented particle method is particularly efficient (see Rem. [3]) . Another com- 
mon choice is v(p) = u max exp (— -j^j , which results in a flux function that has 
one inflection point p* = 2po, and thus can be treated as described in Sect. 16. 11 
For a road network, the road properties (v ma x, Pmax, Po) can differ from 
segment to segment. A proper definition of weak solutions on such networks 
has been presented in [T3] [5] , by introducing "demand" and "supply" conditions 
that have to be satisfied at the nodes. The translation of these conditions to 
the particle method is the subject of current research. 



8.2 Stiff Source Terms 

In Sect. 16 . 21 we have considered the case of simple source terms that are de- 
scribed by the method of characteristics, and that are dominated by the advec- 
tion, so that particle management is well described by the source-free similarity 
interpolation. Here, we outline a case in which the interpolation between parti- 
cles can be used to deal with a dominant source term. 

Many examples in chemical reaction kinetics can be described by convection- 
reaction equations. Such problems are often times governed by a separation of 
time scales. If the advection is scaled to happen on a time scale of 0(1), the 
reactions typically happen on a much faster time scale O(r), where t<1. A 
typical example is the balance law 

«f + (/(«))* - ^(«) , (19) 
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where < u < 1 represents the density of some chemical quantity. The quantity 
is evolved according to a nonlinear advcction term, given by a strictly convex 
flux function /, and a source term. Here we consider a bistable source ip{u) = 
-u(u — l)(it — (3), where < ft < 1. This source term drives the values of u 
towards if smaller than (3, and towards 1 if larger than (3. This example is 
presented for instance in [TT]. Equation (figj) possesses shock solutions as the 
classical Burgers' equation does. In addition, it has traveling wave solutions that 
connect a left state ui«0 with a right state Ur ~ 1- The smooth function that 
connects ul with ur is of width r, and it travels with velocity /'(/3), as shown 
in [5] . The speed of propagation of such solutions is due to the interaction of 
the advection and the reaction term, and many classical numerical approaches 
require the resolution of the wave front in order to obtain the correct propagation 
speed. This can require an unnecessarily fine local mesh resolution of O(r). 

In contrast, the presented particle method possesses key advantages that 
can be used to yield correct propagation speeds and width of the front, without 
actually having to resolve the wave structure. For the bistable reaction kinetics, 
we propose to place and track particles wherever the interpolation crosses or 
touches the zeros of ip (here 0, /3, and 1). In a splitting approach, the particles 
are first moved according to the method of characteristics (fl4|) . Since the source 
is dominant, this step neglects the evolution of the function between particles 
that connect to the value u = (3. Thus, in a correction step, a horizontal motion 
is added to the particles which modifies the area between particles in agreement 
with the source term. The change of area due to the source is calculated by 
using the interpolation (0, and change in area is found by ([6]). While this 
approach is an approximation (the interpolation ([5]) is not the true solution 
between particles), it can be shown that for a Riemann problem between ul = 
and ur = 1, the particle solution converges to a front that moves at speed /'(/3) 
and has a width 0(t). The precise treatment of the transient phase when 
starting with general initial data shall be addressed in future work. 

This approach represents the thin reaction front by exactly three particles. 
Other strategies to treat the stiff source term, which give more control over the 
number of particles on the front, are subject to current research. 

8.3 Operator Source Terms 

Here we consider a conservation law with a source term (I13|) . in which the 
source Su may involve derivatives or integrals of the solution. In these cases, 
the method of characteristics does not apply. Instead, they can be treated by 
fractional steps. In each time step, first the conservation law ut + f(u) x = is 
solved, using the basic particle method described before. Then, in a second step, 
the equation u t = Su is solved. In many cases, the "subgrid" resolution pro- 
vided by the interpolation ([5]) between particles can be employed in this source 
evolution step, either to improve accuracy or to treat the evolution analytically. 
As an example, consider Burgers' equation with a convolution source term 

ut + (^u 2 ) = K * u , (with periodic boundary conditions) (20) 
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which can serve as a model for weakly non-linear waves in a closed, nonuniform, 
slender box [21j . The kernel K embodies the non- uniformity in the spatial do- 
main. For a specific family of kernels K, the mass (J u dx) is always conserved, 
and in addition, the total entropy (here J u 2 dx) is conserved as well, if the 
solution is continuous. The entropy is only diminished at shocks. These prop- 
erties are exactly reproduced by our particle method. In contrast, many other 
numerical methods can diminish the entropy even in the absence of shocks. 

The set of long-term solutions of (f20)) has been conjectured to possess an 
interesting structure. Using our numerical method we are able to examine the 
solution space for evidence of this structure. 

8.4 Higher Space Dimensions 

A scalar conservation law in higher space dimensions 

u t + V • /(«) = , 

with a vector-valued flux function f(u) = (fi(u), . . . , fd(u)), can be approxi- 
mated using dimensional splitting. Each time step consists of the successive 
solution of one-dimensional conservation laws 

«t + (/M)^ = > ( 21 ) 

for i = 1, . . . , d. While this approach is straightforward for Eulerian methods on 
regular grids, Lagrangian approaches have to remesh the solution obtained by 
(p?Tj) . before solving in the direction of the next unknown. In [TS], this approach 
is outlined for front tracking. We propose to use an analogous approach for the 
presented particle method, thus replacing the tracking of shocks by the tracking 
of similarity waves. 

While higher space dimensions can be resolved by dimensional splitting, this 
approach is not fully satisfying, since the required remeshing step spoils one 
of the main benefits of the rarefaction tracking approach: the evolution of an 
exact solution away from shocks. Similarly as for source terms, it may be more 
desirable to use the method of characteristics directly. The characteristic system 

ar = v/(u) 

it = 

yields, as in the one dimensional case, the exact solution on the moving particles. 
However, two fundamental complications arise. First, at a shock, particles may 
not collide, but instead move "around" each other. Second, the definition of 
an interpolation between particles is more difficult. Both aspects are subject to 
current research. 

9 Conclusions and Outlook 

We have presented a numerical method for one-dimensional scalar conservation 
laws that combines the method of characteristics, local similarity solutions, and 
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particle management. At any time, the solution is approximated by rarefac- 
tion and compression waves, which are evolved exactly. Breaking waves, and 
thus shocks, are resolved by local merging of particles. The method is designed 
to conserve area exactly, not increase the total variation of the solution, and 
decrease its entropy. Hence, it is exactly conservative, TVD, and possesses no 
numerical dissipation away from shocks. Shocks can be located accurately by a 
simple postprocessing step. The method is extended to non-convex flux func- 
tions and to simple balance laws. Numerical solutions for convex flux functions 
are second order accurate. 

The approach performs well in various examples. Just a few particles yield 
good approximations with sharp and well-located shocks. This makes the method 
attractive whenever conservation or balance laws arise as subproblems in a large 
computation, and only a few degrees of freedom can be devoted to the numerical 
solution of a single sub-problem. The exact conservation of area, and the exact 
conservation of entropy in the absence of shocks, make the method a natural 
candidate for applications in which conservation of these quantities is crucial. 
We have outlined applications in which these features can make the presented 
particle method advantageous over classical fixed grid approaches. 

Similarly to classical Eulerian and Lagrangian approaches, the presented 
method can be extended to more general equations and higher space dimensions 
using fractional steps. The idea is to treat advection by particles, and approxi- 
mate the solution by similarity waves. Thus, classical ID Riemann solvers can 
be replaced by ID wave solvers. However, fractional step methods incur vari- 
ous disadvantages. A fundamental question is how more general equations and 
higher space dimensions can be treated directly. 
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